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1. Introduction 

The twisted mass formulation of Lattice QCD (see for example references [1, 2]) is being 
studied extensively with Nf = 2 dynamical flavours by the European Twisted Mass (ETM) collab- 
oration [3-7]. In the twisted mass formulation of QCD, the Wilson term is chirally rotated within 
a doublet. For the 2-flavour case, this amounts to a mass degenerate (u,d) quark doublet. In recent 
studies, a partially quenched setup [8] was used to add the strange quark to this formulation. One 
way to include a dynamical strange quark is to add in addition to the strange quark a charm quark in 
a heavier and mass-split doublet. This formulation was introduced in [9] and first explored in [10]. 
We will use that formulation and report on the current status of runs with Nf = 2 + 1 + 1 twisted 
mass quarks conducted by the ETM collaboration. In section 2 we will detail the lattice action and 
the simulation algorithm, while in section 3 we describe our procedure for tuning to maximal twist. 
In section 3.3 we discuss a reweighting procedure allowing to speed up and improve the accuracy 
of tuning to maximal twist. In section 4 we show a set of preliminary results. Finally in section 5 
we describe first results from one level of stout smearing [11]. 

2. Lattice action and simulation algorithm 

In the gauge sector we use the Iwasaki [12] gauge action (with j8 = 1.9 in the ensembles we 
are discussing here). With this gauge action we observe a smooth dependence of phase sensitive 
quantities on the hopping parameter fc around its critical value ((Qnt). The fermionic action for the 
light doublet is given by: 

Si = a 4 [D w [U]+mo,l + mr5^}Xl(x)} (2.1) 

X 

[2], where mo./ is the untwisted bare quark mass tuned to its critical value m cr jt, \l\ is the bare twisted 
quark mass for the light doublet, T3 is the third Pauli matrix acting in flavour space and Dw[£/] is 
the Wilson-Dirac operator. The quark doublet fields that are proportional to the renormalized u and 
d doublet fields in the physical basis are given by: 

yf h » = e te** X l, V? hys = Xie^i^\ (2.2) 

The light doublet twisting angle ft)/ takes the value ft)/ — > ± f as mo./ — m cr i t — > ±0. In the heavy 
sector, as has been shown by Frezzotti and Rossi in [9], a real quark determinant with a mass-split 
doublet can be obtained if the mass splitting is taken to be orthogonal in isospin space to the twist 
direction. The mass term in the heavy sector then becomes: 

S h = a 4 £ {Xh(x) [mo,h + z'Ma75Ti + JU5T3] Xh(x)} . (2.3) 

X 

This time the mass splitting is in the T3 direction and the twist in the Ti direction, which yields a 
similar rotation relation to the physical quark basis: 

Yt" = e-^XH, wt yS = Xue^ . (2.4) 

At maximal twist (|ft)/| = |ft)/,| = |) physical observables are automatically &{a) improved without 
the need to determine any action or operator specific improvement coefficients. The gauge config- 
urations are generated with two different algorithms, depending on whether thin or stout links are 



2 



Status ofETMC simulations with Nf = 2 + 1 + 1 twisted mass fermions 



Siebren Reker 



being used. For thin links a Polynomial Hybrid Monte Carlo (PHMC) updating algorithm [13, 14], 
while for stout links the PHMC setup described in [15] is used. 1 All ensembles described in this 
paper correspond to lattices with spatial extension L/a = 24 and temporal extension T = 2L. Our 
lightest twisted mass quark value (ajj.i = 0.004) corresponds to a pion mass of around 300 MeV 
and a value of m n L around 3.5. 



3. Tuning action parameters 
3.1 Procedure 

Tuning to maximal twist requires to set niQ i and m§ h equal to some proper estimate of the 
critical mass m cr it = >w C rit(j3) [9]. Here we set m Q j = m /, = 1/(2 k) — 4 and, at fixed jj, a and jXg, 
for each value of pii, we tune K* so as to obtain amp C AC,i = (see also eq. (3.1)). Note that the 
untwisted PCAC mass arising from the heavy sector is 0(a) and ft)/, is automatically very close 
to | [10]. Our way of implementing maximal twist can be viewed as the extension to the theory 
with 2 + 1 + 1 dynamical quarks of the optimal critical mass method of Refs. [16]. In this way we 
expect [17] the (absolute) magnitude of lattice artifacts to remain roughly constant as \ii — > 0. The 
PCAC mass in the light sector is defined in the following way: 



mp CA c,i= 2Lx{p a Mpm) ■ (3-D 
where the axial and pseudoscalar bilinears A^j and P/ are 

^ a n,l = XiW^li and P" = Xi75jZt- (3-2) 

The numerical precision at which the condition mpcAC,i — is fulfilled in order to avoid residual 
& (a 2 ) effects when the pion mass is decreased is, for the present range of lattice spacings, | e/jU/ 1 < 
0.1, where e is the small deviation of wipcaci from zero [4, 18]. Note that typically a good statistical 
precision on the estimate of nipcACi is required to satisfy this condition. This computationally 
disadvantageous situation can be remedied through the use of reweighting, as described in section 
3.3. 

The heavy doublet mass parameters and n§ should be adjusted in order to reproduce the 
phenomenological values of the renormalized s and c quark masses. The latter are related to pL c 
and ;Ug via [9] : 

1 Z P 

(» v )r = TjriHo + t^^s) (3.3) 

where the — sign corresponds to the strange and the + sign to the charm. In practice we fix the 
values pL a and jJ.g by requiring the kaon mass to take its physical value and (m c )s. > 10(ra s )R 
(through an estimate of Zp/Zs). 



'Note that in the Monte Carlo simulation we monitor that the minimal eigenvalue of the Dirac operator of the heavy 
doublet always remains positive. 
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3.2 Status 

The tuning status of our ensembles is summarized in table 1 and illustrated in figure 1 . The 
ensembles considered here have the following parameters: j8 = 1.9, a\i c = 0.15 and ajj,§ = 0.19, 
with lattice volumes of 24 3 x 48. 



a\Xi 


0.004 


0.006 


0.008 


0.01 


Kent 


0.163270 


0.163265 


0.163260 


0.163255 



Table 1: Summary of tuning to maximal twist. 
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Figure 1: Tuned runs at yS = 1.9, with a box size of 
24 3 x 48, with a/i a = 0.15 and a/i s = 0.19. From left 
to right the points correspond to a\l\ = 0.004 (red), 
fl/i; = 0.006 (green), aji/ = 0.008 (blue) and afii = 
0.01 (yellow). For all ensembles, »?pcac,i/M/ is within 
the 10% level (see text). 



Figure 2: af n as a function of (am K ) 2 for the four 
ensembles described above. All of these are at /3 = 
1.9, L/a = 24 and T ja = 48. Error bars are plotted, 
but they fall within the size of the dots. 



Tuning to K^it was performed independently for each /I/ value. From table 1 we observe that 
the estimate of K" cr j t depends weakly on Figure 1 shows that for all \ii values the condition 
\e/fXi\ < 0.1 has been fulfilled. 

3.3 Reweighting 

Reweighting towards the chiral limit for (untwisted) Wilson fermions has been recently ex- 
plored in reference [19] (and see also the conference contribution [20] with further references). For 
twisted Wilson fermions one can use reweighting near maximal twist for tuning the last digits of 
the hopping parameter K. Once the update is close to K" cr jt (in the present case for Ak = \ K — Kc r i t \ = 
<^(10~ 5 )) the situation is rather advantageous for an effective reweighting because the gauge field 
configuration is fluctuating between positive and negative values of mpcAC,i- Thus, in order to com- 
plete the tuning one has to just reweight the ensemble to a close-by value of K at which rapcAC,i = 
is realized. 

For calculating the ratio of determinants one could use the methods developed in reference 
[19], but we applied a simple and fast method based on least-square optimized polynomial approx- 
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imations introduced in [21]. In the case of Nf flavours, these polynomials satisfy 

\imP^\x) =x~ Nl/2 , limPi 2 '^) = . 1 x G [s,X] (3.4) 

where the interval [e, X) covers the spectrum of Q*(k)Q(k), with Q the fermion matrix and 
the actual value of n used in the case of the polynomial pW(pW). With these polynomials the 
reweighting factor for fq — > K 2 is 1 / 'det^ij ( K iWi ( k 2)^j»T(Ki)]- A stochastic estimator of the 
reweighting factor can be obtained from a random Gaussian vector as 

expl-T^Q-l)*]}, with ^ = pi 2 2) (jc 1 )Pi 1 1) (K- 2 )P, 1 ( 2 2) (jc 1 ) . (3.5) 

We showed in test runs that with a few stochastic estimators an effective reweighting can be 
achieved for Afc ~ 10~ 5 . We actually applied determinant break-up by a factor of 2, that is, we 
considered two copies of half doublets, both for light and heavy doublets. With this method the 
number of long tuning runs for maximal twist can be reduced. 



4. Results 



In this section we present preliminary results for the light quark mass dependence of the pion 
decay constant and the pion, kaon and nucleon masses (figures 2,3 and 4). The pion masses in 
these ensembles are in the range from 300 to 500 MeV, corresponding to j8 = 1 .9 and a\i\ = 0.004 
to 0.01 (see section 3.2). 
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Figure 3: The kaon (blue) and pion (red) squared 
lattice masses are plotted versus the bare quark mass 
ani. 



Figure 4: The nucleon mass (ootn) as a function of 
(am K ) 2 for the four ensembles described above. A 
subset of the available statistics has been used here. 



We measured ro/a and performed its chiral extrapolation. The results are compared to those 
of our Nf = 2 runs with j3 =3.9 at the same lattice volume in table 2. The vicinity of the two 
values is very encouraging, since it might allow to perform a direct comparison between these two 
ensembles and identify possible effects due to the dynamical strange and charm quarks. 
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•^Gauge 


P 


ro/a in chiral limit 


N f = 2 
Nf = 2 + 1 + 1 


tl-Sym 
Iwasaki 


3.9 
1.9 


5.22(2) 
5.18(5) 



Table 2: Comparison of ro/a in the chiral limit between the Nf = 2 and Nf = 2+ 1 + 1 runs 

5. Stout smearing 

We separately tested the effectiveness and impact on observables of one level of stout smearing 
[1 1] with p = 0.15 (in every direction). In order to perform this comparison, we tuned a run with 
stout smearing at a\l\ = 0.004 to maximal twist. The parameters of this run are compared with 
the one with thin links in table 3. For the simulation with stout smearing the values of the heavy 



p 


traj 


Kcrit 


a\ii 


a\i a 


an s 



0.15 


~5000 
~2500 


0.16327 
0.14552 


0.004 
0.004 


0.15 
0.17 


0.19 
0.185 



Table 3: Comparison of the ensembles with ap/ = 0.004, with (p = 0.15) and without (p = 0) stout smear- 
ing. 

doublet mass parameters p CT and pg were chosen according to the same criterion as described in 
section 3.1. The difference in the values for p CT and pg reported in Table 3 is likely due to the fact 
that the ratio of the renormalization factors Zp/Z$ increases when replacing thin links with stout 
links. Although the statistics are somewhat different between the two runs we can still make a few 
conclusive comparisons. We observe that our estimate of JQnt decreases with stout links. In table 4 
we show a weak dependence of ro/a on a\i[ for simulations with stout links. Finally we observe that 
at fixed bare quark mass a\i\ the pion mass (in physical units) from stout link simulations decreases 
with respect to that from thin link simulations. This finding can be ascribed to an increase of 
the value of the renormalization constant Zp (implying a decrease of the renormalized light quark 
mass). 



a\X x 


0.004 


0.006 


0.008 


ro/a 


4.97(8) 


4.90(8) 


4.99(8) 



Table 4: Comparison of ro/a values for various ensembles with stout smearing at different light twisted 
quark mass value ap/. 



6. Conclusions 

We have presented preliminary results for a set of simple observables obtained from ensembles 
with Nf = 2 + 1 + 1 dynamical flavours of maximally twisted Wilson fermions. We have tuned to 
maximal twist at four values of p/ on 24 3 x 48 lattices at j8 = 1.9. A detailed analysis of the data 
as well as a systematic control of finite size and discretization effects is in progress. 
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